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In this paper we are interested in term structure models for pricing zero coupon bonds 
under rapidly oscillating stochastic volatility. We analyze solutions to the generalized 
Cox-Ingersoll-Ross two factors model describing clustering of interest rate volatilities. The 
main goal is to derive an asymptotic expansion of the bond price with respect to a singular 
parameter representing the fast scale for the stochastic volatility process. We derive the 
second order asymptotic expansion of a solution to the two factors generalized CIR model 
and we show that the first two terms in the expansion are independent of the variable 
representing stochastic volatility. 
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1. INTRODUCTION 

Term structure models describe a functional dependence between the time to matu- 
rity of a discount bond and its present price. Yield of bonds, as a function of matu- 
rity, forms the so-called term structure of interest rates. If we denote by P = P(t, T) 
the price of a bond paying no coupons at time t with maturity at T then the term 
structure of yields R(t, T) is given by 

-R(t,T){T-t) ■ „ D , 

log P(t,T) 



P(t,T) = e- H ^ T ^ 1 - t \ i.e. R(t,T) 



T-t 



(cf. Kwok [7]). Continuous interest rate models are often formulated in terms of 
stochastic differential equations (SDEs) for the instantaneous interest rate (or short 
rate) as well as SDEs for other relevant quantities like e.g. volatility of the short 
rate process. In one-factor models there is a single stochastic differential equation 
for the short rate. The volatility of the short rate process is given in a deterministic 
way. It is assumed to be constant (the Vasicek model or it is a function of the 
short rate itself. In the classical Cox, Ingersoll, and Ross model (CIR) the short rate 
is modelled by a solution to the following stochastic differential equation: 

dr = n(8 — r)dt + a^/rdw, (I) 



2 



B. STEHLIKOVA and D. SEVCOVIC 



where k,9,<j > are parameters representing the rate of reversion, the long term 
interest rate and volatility of the interest rate, respectively (see [5]). By dw we have 
denoted the differential of the Wiener process. Beside these two simple models there 
is a wide range of other term structure models including, in particular, the Chan- 
Karolyi-Longstaff-Sanders model [T], the Hull- White model [5] and many others. 
Based on the assumption made on the form of the short rate process one can derive 
a linear scalar parabolic equation for the bond price as a function of the current 
short rate and time to maturity. 

In the two-factor models there are two sources of uncertainty yielding different 
term structures for the same short rate as they may depend on the value of the 
other factor. Moreover, two-factor models have a richer variety of possible shapes 
of term structures including, in particular, nonmonotone yield curves. The reader is 
referred to the books by Kwok [7] and Brigo and Mercurio [3] for detailed discussion 
on two-factor interest rate modeling. 

There are several ways how to incorporate the second stochastic factor. It is 
reasonable to conjecture that in a financial market the volatility of a fluctuating 
underlying process for the short rate can be fluctuating as well. In the so-called two- 
factor models with a stochastic volatility we allow the volatility to have a stochastic 
behavior driven by another stochastic differential equation. As a consequence of the 
multidimensional Ito's lemma the corresponding equation for the bond price is a 
linear parabolic equation in two space dimensions. These spatial dimensions corre- 
spond to the short rate and volatility. It is well known that the density distribution 
of a stochastic process is a solution to the Focker-Planck partial differential equa- 
tion and can be expressed analytically in the case the volatility undergoes the Bessel 
square root process (see e.g. [S]). The actual value for the stochastic volatility is not 
know in the market. We can just observe its statistical moments like e.g. the mean 
value, volatility, skewness of the volatility etc. Knowing the density distribution of 
the stochastic volatility we are able to perform averaging of the bond price and the 
term structure with respect to the stochastic volatility. Unlike the short rate which 
is known from the market data on daily basis, as it was already mentioned, the 
volatility of the short rate process is unknown. Therefore such a volatility averaging 
is of special importance for practitioners. 

The main goal of the paper is to derive an asymptotic expansion of the bond price 
with respect to a singular parameter representing the fast scale for the stochastic 
volatility process. We derive the second order asymptotic expansion of a solution to 
the two factors generalized CIR model and we show that the first two terms in the 
expansion are independent of the stochastic volatility term. 

The paper is organized as follows. In the next section we present an empirical 
evidence of a short rate process for which the volatility is fluctuating and it has two 
concentration values. Next we discuss a model for statistical distribution captur- 
ing such a volatility clustering. Section 3 is devoted to the asymptotic analysis of 
solutions to the bond pricing equation in the case when the fluctuating volatility 
is rapidly oscillating. We derive explicit formulae for the first three terms in the 
asymptotic expansion of a solution with respect to a small parameter representing 
the fast time scale for rapidly oscillating volatility. 
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Figure 1: Left: Estimates of CIR model's dispersion a 2 from 20-day intervals (3-months 
Treasury bills, 90 intervals starting in January 1990). Right: the density distribution of 
estimates of the dispersion a 2 . 



2. Empirical evidence of existence of volatility clusters and their modeling 

The key feature of the CIR modeling consists in the assumption made on constant 
volatility of the stochastic process (Q]) driving the short rate r. However, in real 
financial markets we can observe a substantial deviation from this assumption. To 
provide an empirical evidence for such a volatility process, we computed maximum 
likelihood estimates of the dispersion for the CIR model for 20-day-long intervals 
using three months treasury bills data. Figured] (left) shows the estimated dispersion 
as a function of time. Higher and lower volatility periods can be distinguished. They 
can be seen also on the kernel density estimates of the values in Figure [T] (right). 

In order to capture such a behavior of the dispersion a 2 we shall consider a model 
in which the limiting density of the dispersion (as t — > oo) has two local maxima. It 
corresponds to the so-called volatility clustering phenomenon where the dispersion 
can be observed in the vicinity of two local maxima of the density distribution (see 
[8]). The desired behavior of the process and its limiting density are show in the 
Figure [2j A natural candidate for such a volatility process is 

dy = a(y)dt + us(y)dw (2) 

having a drift function a(y) such that the differential equation = a(y) has two 
stable stationary solutions. With added stochastic part uj(y)dw of the process, these 
stationary solutions become values, around which the volatility concentrates. Recall 
that the cumulative distribution function G = G(y,t) = Prob(y(t) < y\y(0) = ya) 
of the process y = y(t) satisfying ((2]) and starting almost surely from the initial 
datum yo can be obtained from a solution g — dG/dy to the so-called Focker-Planck 
equation for the density function: 

r)h 1 r) 2 r) 

Tt = 2dy~ 2{w{v)2y9) ~ d~y [a{vY91 ~ 9{V ' 0) = 6{y ~ Vo) (3) 

(cf. Kwok [7]). Here 5(y — yo) denotes the Dirac delta function located at j/o. The 
limiting density g(y) — lim^oo g(y, t) of the process is therefore a stationary solution 



4 



B. STEHLIKOVA and D. SEVCOVIC 




Figure 2: Simulation of a process (left) and its asymptotic distribution (right). 

In [8] one of the authors proposed a model with a property that the limiting den- 
sity is a combination of two Gamma densities. Indeed, let us consider the following 
two mean reverting Bessel square root stochastic processes: 

dyi = K y (9i ~ yi)dt + Vy/y^dwi, i = 1, 2, (5) 

where 9i > 0, 2n y 9i > v 2 > for i = 1, 2, and dwi, dw2 are uncorrelated differentials 
of the Wiener process. Solving the stationary Focker-Planck equation it turns out 
that the limiting distributions of the processes y\ , y 2 are the Gamma distributions 
with shape parameters 2n y /v 2 and 2K y 9i/v 2 . Denote their densities by g\ and g 2 . 

Then gi(y) — Ciy « 2 exp(—-^-y) for y > and gi(y) = otherwise. Here d > 
is a normalization constant such that J R gi(y)dy — 1. Choose a parameter fc € (0, 1). 
Our aim is to construct a process with asymptotic density 

g(y) = kg 1 (y) + (l-k)g 2 (y) : (6) 

corresponding to a convex mixture of densities g\ and g% . In the following theorem 
we see that for the same square root volatility function of the form v^fy it is possible 
to achieve this goal. Drift of the process a(y) can be written as a weighted sum of 
drifts oti(y) — n(9i — y),i = 1, 2, with the weights depending on y. 

Theorem 1 J3 Section 5] Suppose that the drift term a has the form: a(y) = 
w(y)ai(y) + (l-w{y))a 2 (y) where w(y) = kgi(y)/(kgi(y) + (l-k)g 2 (y)) anda.i(y) = 
n(9i — y). Then the stochastic process driven by the SDE: dy = a{y)dt + v^/ydw has 
the limiting distribution g given by the convex combination (0) of densities g± , g 2 ■ 



3. Generalized CIR model with rapidly oscillating stochastic volatility and its 
asymptotic analysis 

The aim of this section is to provide a tool for modeling the effects of rapidly oscil- 
lating stochastic volatility that can be observed in real markets. If the length of the 



On the singular limit of solutions to the CIR interest rate model with stochastic volatiltiy 5 



time scale for dispersion y is denoted by s, the equation for the variable y reads as 

follows: _ 

dy= ^Ml dt+ ^ dWy . (7) 
e \J e 

In what follows we will assume that < e <C 1 is a small singular parameter. Notice 
that the limiting density function g of the stochastic process driven by SDE (J7J) is 
independent of the scaling parameter e > 0. The statement follows directly from the 
stationary Focker-Planck equation ([4]). Concerning structural assumption made on 
the drift function a : R —> R we shall henceforth assume the following hypothesis: 

(A) a is a C 1 function on [0,oo), > ^ limsup < 0. 

V y — >oO y 

Now it is straightforward computation to verify the following auxiliary lemma. 

Lemma 1 Let the drift function a(y) be defined as a mixture of two Gamma limiting 
distributions as in Theorem^ Then the function a satisfies the hypothesis (A) with 



a(0) = k min^!, f? 2 ) and limsup^^ = — n < 



0. 



Next we shall show the limiting density g of the process driven by SDE |(7J) is 
uniquely given by the following lemma: 

Lemma 2 Under the hypothesis (A) made on the drift function a the stationary 
Focker-Planck equation L^g = \-§^s{V9) — ^j( a (y)9) = has a unique solution g 
such that g(0) = for y < 0. It can be explicitly expressed as: 

g(y) = Cy- 1 exp (1 J* ^d^j = Cy^ 1 exp (1 J* 

for y > and g(y) = for y < 0. Here a(y) — (a(y) — a(0))/y and C > is a 
normalization constant such that J °° g(y)dy = 1. 

Proof. It follows by direct verification of the equation. The other linearly indepen- 
dent solution <72 to the equation (J4| has a nontrivial limit 32(0+) > 0. 
In what follows, we denote by a 2 , D > 0, and S the limiting mean value, dispersion 
and skewness of the stochastic process for the y-variable representing stochastic 
dispersion, i.e. 

/>oo />oo -1 />oo 

° 2 =\ yg(y)dy, D= (y-o- 2 ) 2 g( y )dy, S = -y / (y-o- 2 ) 3 g(y)dy . (8) 
Jo Jo Jo 

Notice that D = — f^(£ — °~ 2 )g{£,) d£dy. In the generalized CIR model with a 
stochastic volatility, the instantaneous interest rate (short rate) r will be modelled 
by the mean reverting process of the form |T]) where the volatility of is replaced by 
a square root of a stochastic dispersion y, i.e. 



dr = k(9 — r)dt + ^Jy\frdw r . 



(9) 
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The differentials of the Wiener processes dw y and dw r are assumed to be uncorre- 
cted throughout the paper, i.e. E(dw y dw r ) = 0. Then the corresponding partial 
differential equation for the bond price P £ — P £ {t,r,y) has the following form: 



dP £ , ,„ , r , , i ^,dP £ 1 d 2 P e 



+ ( K (9-r)-\ 1 (y,r)r2^)— + -ry— 1 -~rP £ (10) 



dt y y ' lvy ' ' yi " dr 2 u dr 2 

1 / - , , ^dP £ \ 1 / , ,dP £ v 2 yd 2 P £ 
-M{y,r)v^y——)+-[a(y)- 



yfe \ dy J e \ dy 2 dy' 

(t,r,y) E Qt = (0,T) x R + x R + , where Ai,A 2 are the so-called market prices of 
risk (cf. [7J Chapter 7]). By a solution P £ to (1 1 1 1) we mean a bounded function 
P £ E C 1,2 (Qt) H C(Qt) satisfying equation (jlip on Qt- Concerning the structural 
form of market prices of risk functions Ai, A2 we shall suppose that 

Ai(t,r,y) = XiVrVv, M*, r, y) = A 2 ^y 

where Ai,A2 € -R are constants. It is worthwilc noting that the latter assumption 
is not restrictive as the original one-factor CIR model assumes such a form of the 
market price of risk (cf. Kwok 7J). We shall rewrite PDE (TTT]) in the operator form: 

(e- 1 £ + e- 1/2 £i+£ 2 )P £ = 0, (11) 
where the linear differential operators £q,£i,£2 are defined as follows: 



d v 2 y d 2 d 

C ° = a{v) d-y + —W 1= 2VV d-y 

d d Id 2 



£2 = ^ + (^-^)- A i^)^: + 2 r ^- r - ( 12 ) 

Next we expand the solution P e into Taylor power series: 

00 

P £ (t,r,y)=J2^Pj(t,r,y) (13) 
j=o 

with the terminal conditions Po(T,r,y) = l,Pj(T,r,y) = for j > 1 at expiry 
t = T . The main goal of this paper is to examine the singular limiting behavior of 
a solution P £ as e — ► + . More precisely, we shall determine the first three terms 
Pa, Pi, P2 of the asymptotic expansion (IT5]) . We shall henceforth denote by (ip) 
the averaged value of the function tp E C([0, 00)) with respect to the density g, i.e. 
(if)) = J °° ip(y)g(y) dy. We shall also use the notation (£ 2 ) standing for the averaged 
linear operator £ 2 , i.e. 

<£ 2 ) = | + - r) - A 1W 2 ) |- + ia 2 r^ - r (14) 



Lemma 3 Let ip E C 1 ([0, 00)) be such that £gip is bounded. Then (£o^P) = 0. 
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Proof. Notice that the operator £q is the adjoint operator to the linear operator 
Cq with respect to the L 2 -inner product {tp,4>) = J i } {y) ( P{y) dy. It means that 
(Coip) = (£oip,g) — (-0, £oS) = because the density g is a solution to Eq. (01. (} 
The following lemma will be useful when computing higher order term in series 
expansion (fT5| . 

Lemma 4 Let F € C([0, oo)) be such that (F) = 0. Then, up to an additive 
constant, there exists a unique solution tp £ C 2 ((0,oo)) n C([0, oo)) to the non- 
homogeneous equation C,Qip = \F. Its derivative ^ is given by 

dtp 1 f v 

oy yg{y) J 

Moreover, {Citp) — X2V J °° F(y)yg(y)dy. In particular, if ip is a solution to the 
equation Cq^P — then tp is a constant function with respect to the y-variable. 

Proof. Using equation for the limiting density g and inserting ^ into the 

2 

operator £0 we obtain that ip is a solution to the equation C^ip = ^-F. Other 
independent solutions are not continuous at y = 0. The formula for (jCiip) follows 
from the definition of the operator L\ by applying integration by parts formula. <0 

Now we proceed with collecting the terms of the power series expansion of (jlljl . 

• In the order e _1 we have LqPq = 0. According to Lemma[4]we have Pq = Po(t, r), 
i.e. Pq is independent of the y-variable. 

• In the order e -1 / 2 we have LqP\ + C\Pq = 0. Since Pq = PQ(t,r) we deduce 
CiPq — and so CqP\ — 0. By LemmalU Pi = P L {t,r) is independent of y. 

• In the order e° we have £0^2 + £>\Pi + £2^0 = 0. Since Pi — Pi(t,r) we have 
C\Pq = 0. Hence £0-^2 + C2P0 = 0. Taking the average (.) of both sides of the 
latter equation we obtain (£0^2) + (£2^0) = 0. By Lemma [3] and the fact that Pq 
is independent of y-variable we conclude (£2)^0 = (£2^0) = 0. Therefore Pq is a 
solution to the classical one-factor PDE equation for the CIR model satisfying the 
terminal condition Pq{T,t) = 1 for any r > 0. It is well known that the solution 
Pq = Po(t, r) to the equation (£2)^0 — is given by the explicit formula: 

PQ(t,r)=A (t)e- B ^ r , (15) 

where A' = k9B and £' = («+ \io 2 )B + ^-B 2 - 1, A (T) = 1, BIT) = 0, i.e. 

_ / 20e (^)(T-t)/ 2 x y _ 2(e*Cr-q - 1) 

0U U0 + 0)(e^ T -*)-l) + 2^ ' w (0 + V)(e^ T "*) -l) + 20' 

= k + Aicr 2 , = V^ 2 + 2ct 2 (cf. Kwok [II Chapter 7]). Since (£ 2 )P = we have 
-£ 2 P = ((£2) - £2) Po = (o- 2 - y)f{t)re- B ^ r where 

f{t) = {\iB{t)+ l -B{t) 2 )AQ{t). 
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Hence L0P2 = —C2P0 = [p 2 — y)f(t)re B w r . According to LemmaHwe have 
ftp o 1 f y 

W = - V ,mre-^>H{y), H(y) = ^ (f - a*)^. (16) 

• In the order e 1 / 2 we have £op3 + AP2 + £2 -Pi = 0. Since (£op3) = we have 
(C1P2} + (£ 2 Pi) = 0. The function Pi = P\{t,r) is independent of the y-variable 
and therefore (£2) -Pi = (£2 Pi) = — (£ip2)- By Lemma|4]we have 

£iP 2 - -^f(t)re- B ^ r yH(y), -(£iP 2 ) = K.f^re^^ , 

where ^ = Jq(£ ~ v 2 )g(0d£dy = ^r D is a constant (see ©). Notice 

that the constant K\ and the function /(i) depend on the first two moments a 2 and 
D of the stochastic dispersion only. Equation (£ 2 )Pi = (£2 Pi) = — (£iP 2 ) reads as: 

^ + («(* - r) - W) ^ + - rP! = ^/(t)re-^. (17) 

The solution Pi satisfying the terminal condition Pi (T, r) = for r > can be 
found in the closed form: 

Pi(t, r) = (A w (t) + A 11 (t)r)e- B ^ t > (18) 

where the functions Aio(i), Au(t) are solutions to the system of linear ODEs: 

A' n (t) = (K9B{t) + K + \ lC r 2 +<r 2 B(t)) A lx (t) + K x f{t), (19) 
A' w (t) = K0B(t)A w (t)- K0A u (t), 

with terminal conditions Aio(T) = 0, An(T) = 0. We can analytically and also 
numerically compute A%o, An in a fast and accurate manner. This way we have 
obtained the term P\(t, r). In Figure [3] examples of numerical approximation of the 
term structure R £ (T — t, r) = — - log(P e (P — r, r, .)) corresponding to the second or- 
der expansion of the averaged value of (P e (i,r, .)), P £ (t,r, y) w Po(i, r) + yfeP\(t, r). 
We plot term structures starting from the short rate r = 0.03 (left) and r = 0.031 
(right) for parameters k = 5, 6 = 0.03, n y = 100, u = 1.1832, d 1 = 0.025 , 2 = 0.1, 
jfc = 1/3, Ai = -1, A 2 = -100 and e = 0,0.001,0.01 (black, red and blue curves). 

Having Pi and we can compute the term £iP2 + £ 2 Pi- With regard to 
Lemma |4] equation £oPj = —C1P2 — C2P1 then yields a formula for and 

(£iP 3 ) = + ^)/(*)re-^ + P(-Ai,^ + I 

where the constant K3 — J °° £ 3 H(£)g(£)d£ = —^SD^ — a 2 D depends on the first 
three statistical moments of the stochastic dispersion. 

• In the order e l we have £oPi + £ip3 + £ 2 P2 = 0. Proceeding similarly as before 
we have (£oPi) = and therefore 

(£iP 3 ) + (£ 2 P 2 ) = 0. (20) 
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Figure 3: The approximate term structure R e = R e (T — r, r) based on the first two 
leading terms of the bond price P E w Po(T — r,r) + y/ePi(T — T,r) starting from the 
short rate r — 0.03 (left) and r = 0.031 (right) for several values of the singular parameter 
e = 0, 0.001, 0.01 (black, red and blue curves), resp. 



We decompose the function P 2 {t, r, y) in the form 

P 2 (t,r,y)=P 2 (t,r)+P 2 (t,r,y), (21) 

where P 2 is the averaged value of P 2 and P 2 is a zero mean fluctuation, i.e. (P 2 ) = 0. 
As P 2 does not depend on y, we have = Taking into account (C 2 P 2 ) = 
we obtain 

P 2 (t,r,y) = -^f(t)re- B ^ (J^H(0d^-K 2 \ 

where K 2 — g(s) H^d^ds is a constant and the function H is given by 
(|16|l . Now we can use decomposition ([21]) to evaluate (C 2 P 2 ). We have (£ 2 P 2 ) = 
(C 2 (P 2 + P 2 )) = (£ 2 )P 2 + {C 2 P 2 ) because P 2 is independent of y. Next we can 
determine (C 2 P 2 ) in the following form: 

where K4 = J °° H {^)d^(y — a 2 )g{y)dy . It is worthwile noting that both constants 
K 2 ,K^ depend on all nontrivial statistical moments of the stochastic dispersion. 
Equation (|2"0|) then becomes 

(C 2 )P 2 = -(C 2 P 2 ) - (AP 3 ) = (a(i) +b(t)r + c(t)r 2 )e- B ^ r , P 2 (T,r) = 0, 

which is a partial differential equation for P 2 = P 2 (t,r,y) with a right hand side 
which can be explicitly computed from already obtained results in the closed form: 

P 2 (t,r) = (A 20 (t) + A 21 (t)r + A 22 (t)r 2 )e~ B ^ r (22) 

where the functions A 2 o, A 2 i,A 22 are solutions to a linear system of ODEs. We omit 
details here. 

In summary we have shown the following main result of this paper: 



10 



B. STEHLIKOVA and D. SEVCOVIC 



Theorem 2 The solution P £ = P e {t, r, y) of the generalized CIR bond pricing equa- 
tion ill]) with rapidly oscillating dispersion can be approximated, for small values of 
the singular parameter < £ « 1, by P E (t, r, y) ~ Po{t, r) + y/sPi (t, r)+eP2(t, r, y) + 
0( £ l). 

The first two terms Pq,Pi are independent of the y-variable representing unob- 
served stochastic volatility. They depend only the first two statistical moments (mean 
value and dispersion) of the stochastic dispersion and other model parameters. 

The next term in the expansion P2 nontrivially depends on the y-variable. P2 as 
well as its averaged value (P2) depends also on all nontrivial statistical moments of 
the stochastic dispersion. 

The terms P$,P\,P2 can be evaluated by closed-form formulae H5\) . H8\) . (2ty). 
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